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Convergence  of  a  Two-Stage  Richardson  Iterative  Procedure 
for  Solving  Systems  of  Linear  Equations. 

*  t 

Gene  H.  Golub   and  Michael  L.  Overton 


0 .  Introduction 

Consider  the  problem  of  solving  a  system  of  linear  equations 

Ax  =  b  (0.1) 

by  an  iterative  method,  i.e.  generating  a  sequence  of  approximations 
{x  }  such  that  lim  x   =  x.   Frequently  it  is  useful  to  introduce  a 
splitting 

A  =  M  -  N 

where  systems  of  the  form 

My  =  c  (0.2) 

may  be  solved  much  more  easily  than  the  original  system  (0.1)  .   When 
(0.1)  arises  from  the  discretization  of  a  partial  differential  equation, 
such  a  splitting  is  often  natural,  with  M  and  N  corresponding  to  sepa- 
rate terms  in  the  differential  equation.   The  iterative  method  used   to 
solve  (0.1)  can  then  be  "preconditioned"  by  M,  i.e.  designed  so  that 
each  step  of  this  "outer"  iteration  involves  solving  a  system  of  the 
form  (0.2).   Sometimes  it  is  desirable  to  solve  these  systems  (0.2)  by 
"inner"  iterative  procedures.   In  this  paper  we  consider  using  the 
second-order  Richardson  method  for  the  outer  iteration,  and  show  how 
the  rate  of  convergence  of  this  iteration  depends  on  the  accuracy  re- 
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quired  for  the  inner  iterations.   To  our  knowledge  analysis  of  this 
particular  method  has  not  been  attempted  previously.   See  Nichols  (1973) 
for  an  analysis  of  more  general  schemes  using  inner  and  outer  iterations 
to  solve  linear  systems.   Other  related  papers  on  the  solution  of  linear 
or  nonlinear  problems  by  two-stage  methods  include  Gunn  (1964),Nicolaides 
(1975)  ,  Pereyra  (1967)  and  Dembo,  Eisenstat  and  Steihang  (1980)  .   See 
Golub  (1962)  for  a  study  of  round-off  error  for  the  Richardson  method, 
which  we  do  not  explicitly  consider  here. 

One  motivation  for  our  work  is  that  the  two-stage  method  can  be 
effectively  used  to  solve  nonsymmetric  systems  where  M  is  s\Tmnetric 
positive  definite  and  N  is  skew-symmetric.   In  this  case  the  symmetric/ 
skew- symmetric  splitting  A  =  M-N  is  used  to  precondition  the  outer  itera- 
tion, and  a  symmetric  splitting-  M  =  M, -M_  can  be  used  to  precondition 
each  inner  iteration,  using  a  direct  method  to  solve  a  system  of  the 
form  My  =  c  at  each   step  of  each  inner  iteration.   (See  Manteuffel 
(1977)  for  alternative  approaches  for  nonsymmetric  systems.)   Numerical 
results  from  applying  the  Richardson  method  to  such  nonsymmetric  systems 
are  given  in  Section  3.   We  also  present  numerical  results  using  the 
conjugate  gradient  method  for  the  outer  iteration.   Analysis  of  the 
latter  procedure  would  be  very  interesting,  but  this  seems  more  diffi- 
cult. 

We  use  II  •  II  to  denote  the  Euclidean  vector  and  matrix  norm,  defined 


by 


II  /    T    ^  1/2  ii^ii  llBxl 

xli    =    (x   x)    '^  ,       llBll    =  max  j, — jr- 

11x11*0   '^^" 


1.   Symmetric  positive  definite  systems. 

Let  us  first  assume  that  both  A  and  M  are  sv-mmetric  and  positive 
definite.   Consider  the  following  iterative  method. 

Method  1. 
Choose  positive  scalar  parameters  5  ,  a  and  lo,  with  0  <  5  <  1. 
Choose  initial  vectors  x   and  x,.   For  k  =  1,2,...,  define 

^k+1  =  ^k-1  ^  '"^^^k^^k-^k-l^  (1-^^ 


where 


MZj^  =  ^k  +   qj,    '  (1-2) 

r^  =  b  -  AXj^  , 
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and 

llqj^ll     <ollrj^ll.  (1.3) 

If  6  =  0  then  Method  1  reduces  to  the  second-order  Richardson 
method.   The  conditions  on  the  linear  operator  and  on  the  parameters 
a  and  co  under  which  convergence  is  guaranteed  are  well  known  in  this 
case  (see  Golub  and  Varga  (1961)).  _  When  6  >  0  the  meaning  of  (1.2)  and 
(1.3)  is  that  the  system  .  ..-^^,.,,^^, 

Mz  =  r^^  --.-    :i---      --:-  (1.4) 

is  being  "solved"  by  an  inner  iterative  procedure,  which  is  terminated 
when  the  residual  norm  " q,  "  has  been  sufficiently  reduced.   If  the 
outer  (main)  iteration  is  converging  to  the  solution  of  (0.1),  the 
associated  residuals  ^^r  }  are  converging  to  zero,  and  hence  z  =  0  is  a 
reasonable  starting  point  for  each  inner  iteration.   Equation  (1.3) 
then  specifies  that  each  inner  iteration  must  reduce  its  associated 
residual  norm  by  a  factor  of  6 .   Note  that  the  early  inner  iterations 
solve  (1.4)  with  low  absolute  accuracy,  while  the  later  systems  are 
solved  with  high  accuracy.   There  is  no  restriction  on  the  type  of 
method  used  for  the  inner  iterations. 


1 . 1   Convergence  analysis 

Let  us  now  analyse  the  convergence  of  Method  1.   Let  the  error  at 
each  step  of  the  outer  iteration  be 


e,    X    X.  , 

k         k 

where  x  is  the  solution  of  (0.1) .   Note  that  r,  =  Ae,  for  all  k.   For 
k  =  1,2,...,  we  have : 

^k+1  =   %-l    -    '"^"^k  ^  ^k-1  -  ^k^ 

=  ojKe^  +  (l-(^)e^^_-^  +  pj^  (1.5) 


where 


and 


p^   =    -waM   q^ 


K  =  I  -  aM  "''A.  (1.6: 
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Now  K  is  similar  to  a  symmetric  matrix,  and  we  can  define  its  eigenvector 
decomposition  by 


K  =  M-^/2(I-aM-l/2^M-l/2^Ml/2  =  ^-l/2^,^V/2 


:i.7) 


where  Z  is  a  diagonal  matrix  of  eigenvalues  {a.},  j  =  1, 
orthogonal.   Let 


,n,  and  V  is 


„T„l/2      ,  -  ,rT,,l/2         „T  -1/2 

e,  =  V  M    e,  and  p   =  V  M    p   =  -(i)aV  M     q   . 


(1.8) 


Then  from  (1.5)  we  obtain  the  diagonalized  system  of  difference  equations 


^k+1  "  "^^^k  "^  ^-^""^^^k-l  "^  ^k  '  '^  ""  1/2,. 


(1.9) 


At  this  point  we  state  a  lemma  which  can  be  proved  using  the  standard 
theory  of  difference  equations. 

Lemma  1 .   Consider  the  inhomogeneous  difference  equation  in  E,    : 


?k+l  =  e^k  ^  ^^k-1  -^  ^k'  '^  =  ^'2' 


(i.io; 


The  solution  to  (l.lO)  is  given  by 


k-1 


^k  =  ^k  ^1  ^  ^\-l  ^0  ^   J^  \-£^iL 


where 


^1  -  ^2 
H  -  ^2 


kA 


k-1 


if  X^   *  ^2 


if  A,  =  X2 


and  where  A,  and  A-  are  the  roots  of  the  characteristic  polynomical 


A-eA-Y=0 


(Note  that  if  the  coefficients  of  (l.lO)  are  real,  A-,  and  A-,  may 
complex  but  6   is  real) .  Q 

It  follows  from  Lemma  1  that  the  solution  to  (1.7)  is 


be 


k-1 


^k  =  \^1  ^  (l-^)Vl^O  ^   Ji  2^-A'   ^  =  1/2,...,      (.1.11) 
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where  S   is  a  diagonal  matrix  with  j    diagonal  element 


^^k^:, 


kA 


k-1 


l,j 


if  A,  .     *    X^     . 
1,3  2,3 


if  X,  .  =  X_  . 
1/D     2, J 


(1.12; 


and  where  X^     .    and  X„  .  are  the  roots  of 
1/D       2,3 


A   -  Loa .  A  +  (lo-1)  =0,    j  =  l,...,n. 


(1.13) 


Now  let 


p  =  max    (max  (  |  A,   .  |  ,  |  A_  . 


(1.14) 


It  can  easily  be  shown  (by  induction  on  k)  that 


S,  II  E  max   1  (S,  )  .  .   <  kp 

^      i<j<n    k  ::'  - 


k-l 


No 


te  also  from  (1.13)  that  the  product  of  the  roots  A-,   .  A^  •  is  equal 


to  co-1,  so  p   >_  |lo-1|.   We  therefore  have  from  (1.11)  that 

k-l 


e.ll  <  kp^"^lle,  II  +  (k-l)  p*^ lie -II  +    I       (k-£)  p 
^      ~  ^  "      5,=1 


k-Ji-1 


P^ll  ,  k  =  1,2, 


At  this  point  we  need  a  lemma  relating  Hp^H  to  He,  II 


Lemma  2.   Assume  that  (1.3)  holds.   Then 


Pk"  -  ^"^k"  '    k  =  1,2, 


(1.15; 


where  e  =  6  wall  m'-^'^^II  II  AM""'"'^^I 


(1.16) 


Proof .  By  (1.8)  and  (1.3)  we  have; 


p^ll    <    ojallM    -^"^^ll  II  q^ 


<   6ajallM    ^^^H  H  r. 
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-1/2  -1/2      T    1/2 

<   6coallM    ^     II  HAM      '^  W   M^     e. 


from  which   the    result    follows.     []] 

Continuing   with    the    error    analysis,    we    can    substitute    e  II  e    II    for 
llpjl     in    (1.15),    and    rewrite    it    to    obtain 

k-1 

He,  II    -    £       y       (k-Ji)p^"^"^    lleJI  <      kp'^""'-    He,  II    +    (k-l)  p'^    lle^ll, 
£=1  ...:,,..-,..". 

k    =    1,2, ...,m,     (1.17) 

where  m  is  introduced  to  indicate  the  last  computed  iterate  x  .   This 

m 

system  of  m  linear  inequalities  can  be  written  using  matrix  notation  as 


(I  -  eL) 


m 


<  s 


where  the  non-negative  strictly  lower  triangular  matrix  L  and  the  vector 


m 


s  are  defined  accordingly.   Now  L   =  0  so 


/-r      T\~l         -r    _L      T    j_      2^  2  m-l,.m-l 

(I-eL)         =I+eL+eL      +...+e         L 


Thus  (I-eL)   _  is  also  a  non-negative  matrix  and  hence 


m 


<  (I-cL)  "^s 


(1.18) 


T  -1 

Let  us  define  t  =  [t,,...,t  1   =  (I-eL)   s.   By  definition,  t,  satisfies 

1       m  k 


k-1 
T,  =  kp^~^  lleJI  +  (k-Dp^lle^ll  +  e   T   (k-il)  p'^"^"^ 


^l' 


k  =  1,2, ...,m. 


The  first  equation  defines  t  ,  =  II  e,  II  .   For  convenience,  we  now  define 

T„  =  -lle^ll  ,  and  it  then  follows  from  Lemma  1  that  t,  is  the  solution 
0       0  k 

to  a  difference  equation  of  the  form  (l.lO)  with  inhomogeneous  term 
ET   and  characteristic  polynomial  with  a  double  root,  i.e. 
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This, however ,  can  be  viewed  as  the  homogeneous  equation 

and  hence,  using  Lemma  1  again,  has  solution 
k     k  k-1     k-1 

T,   = T,   -   P   .   T-   ,  1.19 

k    v^  -  v^   1    ^     ^1  -^2-  '    ° 


where  v,  and  v_  are  the  roots  of 

X^  -  (2p+e)  A  +  p^  =  0. 
If  we  define  p  =  max(v,,v„)  we  have 

p  =  P  +  f  +  (P£  +  ■^)  ^""^  (1-20) 

=  p  4-  pl/2,1/2  ^  1^  o(e'^^),    if  e    «    1.    , 

Using  (1.18)  and  (1.19)  we  now  have  a  condition  which  guarantees  the 
convergence  of  the  error  norms  { II  e,  II  }  to  zero,  namely  p  <  1.   More 
specifically,  we  have  the  following  result: 

th  j- 
Theorem  1.  The  error  norm  lie,  II  associated  with  the  k  iterate  of 
k 

Method  1  is  bounded  by 

lle^^ll  <  kp^"^  lli^ll  +  (k-l)p'^  lle^ll  , 

where  p  =  p  (6  ,a,La)  is  given  by  (1 .  20)  ,  ( 1 .  16  )  ,  (1 .  14  )  ,  (1 .  13)  and  (1.7)  .  Q 
Note  that  if  6  =  e  =  0 ,  we  have  p  =  p  and  Theorem  1  reduces  to  the 
standard  convergence  property  of  the  second-order  Richardson  method 
(see  Golub  and  Varga  (1961) ) .   We  also  note  that  it  is  possible  to 
choose  the  first  iterate  x,  in  such  a  way  that  the  error  bound  of 
Theorem  1  is  somewhat  reduced. 
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1 . 2   Optimal  parameters. 

There  are  three  free  parameters  in  Method  1,  namely  5,  a  and  lo. 
Let  us  first  suppose  that  6  =0,  so  that   the  inner  systems  (1.4)  are 
solved  exactly.   In  this  case  the  optimal  choices  of  a  and  to/  i.e.  those 
that  minimize  p  in  (1.14),  may  be  described  using  Young's  theory  of 
successive  overrelaxation.   Let  ;  - 

I  0^  I  =  max    ]  a.  |  j  -       -- 

where  {a.}  are  the  eigenvalues  of  K  =  I-oM   A,  as  before.   Given  any  a» 

-'  1 

the  optimal  choice  of  lo  is  obtained  by  choosing  A,    =  X„     =  ^coo  „  in 

(1.13) ,  i.e. 

oj  =  ^        .  (1.21) 

1  +  A-a^ 

It  can  be  verified  that  this  choice  results  in 

P  =  1^1,  J  =  1^1, jl  =  l^2..jl  =  I'^l'^il  =  "^^^ 

for  all  j,  j  =  l,...,n.   Note  that  (1.21)  implies  1  <  oj  <  2 .   See 
Varga  (1962, p. 109)  or  Young  (1971, p. 171)  for  the  proof  that  (1.21)  is 
optimal,  and  see  Golub  and  Varga  (1961)  for  a  description  of  the  connec- 
tion between  the  Richardson  and  successive  overrelaxation  methods. 

Since  p  =  -=-io|aj|,  and  lo  satisfies  (1.21),  it  is  clear  that  a 
should  be  chosen  to  minimize  | a „ | ,  the  spectral  radius  of  K.   Let 

y     and  y  .   be  respectively  the  maximum  and  minimum  eigenvalues  of 
max       mm        r-        j  n 

M   A.   Clearly  \o^\     is  minim.ized  by  specifying 

1  -  ay  .   =  -  (1-ay    )  , 
mm         max 


I.e. 


c  =  ; (1  22) 

y     +  y  .  V  J- .  ^z; 

max     mm 


and 


y    -  y  • 

I  „  I     max    mm   ^  , 

I^J  =  ^ +  y  ■    ^  ^  • 

max    ^mm 

When  we  permit  6  >  0  the  optimal  choice  of  parameters  is  more  com- 
plicated.  Our  aim  is  to  minimize  the  amount  of  work  needed  to  obtain 
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a  solution  of  prescribed  accuracy.   A  reasonable  measure  of  the  total 
amount  of  work  required  for  generating  x, ,...,x   is 

m 

W  =   I    m 
k=l    ^ 

where  m,  is  the  number  of  iterates  required  in  the  k    inner  iteration 
k 

to  achieve  (1.2)  and  (1.3).   This  measure  is  particularly  appropriate 
when  each  step  of  every  inner  iteration  involves  solving  another  system 
of  equations  by  a  direct  method.   Now  let  y  tie  a  measure  of  relative 
accuracy  prescribed  for  the  outer  iteration,  e.g.  specifying  II  e  II  < 
Y lie  II.   Roughly  speaking,  Theorem  1  says  that  the  error  is  reduced  by 
approximately  a  factor  of  p  at  each   step  of  the  outer  iteration  (this 
neglects  II  e,  II  and  a  factor  of  (k+l)/k).   With  this  assumption  the  number 
of  iterates  required  for  the  outer  iteration  is  approximately 


-log  Y 
m  " — '- 


-log  p  (6  ,  a  ,(jo) 


Now  let  us  suppose  that  the  iterative  method  used  for  the  inner  itera- 
tion is  such  that  the  associated  residual  is  reduced  by  about  a  factor 
of  e  at  each  step.   Assuming  a  starting  point  z  =  0  for  each  inner 
iteration  we  have 

"log  5        ,    T      _ 
m.  ~    — T — - — ^     ,   k  =  1 ,  .  .  .  ,  m . 
k    -log  e 

Under  these  assumptions  the  total  amount  of  work  W  is  about 

-log  Y  ,  -log  5 


VJ 


-log  9    -log  p  (6  ,a,Lo) 


Now  Y  and  9  are  fixed  so  a  reasonable  goal  in  choosing  the  parameters 
5  ,  a  and  u)  is  to  minimize 


-log  p  {6   ,a  ,aj) 


Naturally  this  is  difficult  to  do,  but  the  formula  gives  some  insight. 
As  6  ->■    0,  W  ^  °°,  indicating  that  solving  the  inner  iterations  too 
accurately  is  very  expensive.   On  the  other  hand,  if  6  is  allowed  large 
enough  so  that  p  ->  1 ,  we  have  W  -^  =° ,  indicating  that  the  outer  iteration 
is  not  convergent.   The  optimal  value  of  6  ,  given  a    and  lo  ,  is  somewhere 
in  between  these  extremes. 

The  optimal  values  of  a  and  co  are  equally  complicated,  since  p  de- 
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pends  on  them  through  both  p  and  e,  but  it  seems  likely  that  (1.21)  and 
(1.22)  would  be  good  choices.   One  reasonable  way  to  choose  the  para- 
meters in  practice  is  to  make  several  numerical  experiments.   Fortunate- 
ly, once  an  adequate  choice  of  parameters  is  determined  for  a  particular 
problem,  whole  classes  of  related  problems  can  usually  be  solved  effi- 
ciently with  the  same  parameter  values.   Another  interesting  and  promi- 
sing idea  is  to  try  to  choose  6  dynamically,  as  the  outer  iteration 
proceeds . 


2 .   Symmetric/skew-symmetric  splittings. 

Method  1  can  be  used  to  solve  (0.1)  if  A  is  nonsymmetric  but  with 
a  positive  definite  symmetric  part.   We  define  the  splitting  A  =  M-N 

by 

IT  IT 

M  =  f (A  +A) ,  N  =  -(A  -A) , 

where  M  and  N  are  respectively  called  the  symmetric  and  skew-symmetric 
parts  of  A,  and  M  is  positive  definite.   The  matrix  K  is  given  by 

K  =  I-aM~  A  =  (l-a)I  +  aM"  N. 

Now  M~  N  is  similar  to  the  skew-symmetric  matrix  M~  ^  NM     ,  and  thus 
its  eigenvalues,  say  [Ik.},    j  =  l,...,n,  are  imaginary.   The  eigen- 
values of  K  are  thus  given  by 

a  .  =  (1-a)  +  ^a  K  .  .  (2.1) 

Let  us  now  investigate  the   choice  of  optimal  parameters  when 
6=0.   Young  (1971, p. 196)  gives  formulae  for  the  optimal  choice  of  tu 
and  the  resulting  value  of  p  in  the  case  when  the  eigenvalues  of  K  are 
complex  but  confined  to  the  rectangle  defined  by  the  real  and  imaginary 
parts  of  the  largest  eigenvalue  in  modulus,  a  property  which  (2.1) 
satisfies.   It  can  be  verified  that  p  is  minimized  when  a    in  (2.1)  is 
chosen  by 

a  =  1.  (2.2) 

With  this  choice  the  eigenvalues  [a.]    are  imaginary,  and  the  optimal 
choice  of  to  reduces  to  (1.21).   Notice  that  now  X^    ^    =    X^    ^    is    imagi- 
nary and  0  <  oj  <_  1 . 
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Having  motivated  the  choice  (2.2) ,  we  now   permit  6  >  0  and  con- 
sider the  convergence  of  Method  1  using  the  symmetric/skew-symmetric 

splitting,  with  a  =  1.   All  of  the  error  analysis,  including  Theorem  1, 

T       H 
applies  to  this  case,  provided  only  that  we  replace  V   by  V   in  (1.7) , 

where  V  is  now  a  complex  unitary  matrix  (and  Z  is  imaginary) . 


3 .   Numerical  Experiments 

Consider  the  differential  equation: 

-Au  +  (au)   +  au   +  (bu)   +  bu"  " +"  ciT  =  f  (3.1) 

XX       y     y 

where  u,a,b  and  f  are  functions  defined  on  the  unit  square:  0  <  x  <  1, 


0  <  y  <  1.   Approximating  (3.1)  by  second-order  finite  differences  on 
a  grid  with  t  interior  poin 
linear  system  to  be  solved: 


2    2 
a  grid  with  t  interior  points  in  each  direction  results  in  an  t  x  t 


(M-N)u^  =  d.  (3.2) 

T 
Here  u   is  the  solution  to  the  discretized  problem,  N  =  -N   xs  a  skew- 
symmetric  matrix  corresponding  to  the  first-order  terms  in  (3.1) , 
M  =  M   =  M, -M   is  a  symmetric  matrix,  with  M   the  negative  discrete 
Laplacian  operator  and  M   a  diagonal  matrix  corresponding  to  the  zero- 
order  term,  and  d  corresponds  to  f,  modified  to  account  for  the  boundary 

conditions.   We  consider  the  particular  differential  equation  whose 

2   2 
solution  is  u(x,y)  =  exp(x  +y  ),  with  coefficients  a(x,y)  =  b(x,y)  = 

2   2  2    2 

5  exp(x  +y  ),  c(x,y)  =  10  exp(3x  +3y  ),  and  with  nonhomogeneous  Dirichlet 

conditions  imposed  on  the  boundary  of  the  square. 

Table  1  gives  the  result  of  solving  (3.2)  by  Method  1.   Note  that 

the  symmetric/ skew- symmetric  splitting  was  used  for  the  outer  iteration, 

as  discussed  in  Section  2.   The  parameters  were  a  =  1  and  co  given  by 

(1.21)  (see  below).   The  initial  vector  x   was  set  to  zero,  and  x   was 

obtained  by:  x,  =  z  +x  ;  Mz   =  r  +q  ;  llq-ll  £  5  H  r  11  .   The  termination 

criterion  was  II  r  11  <  10~  .   Each  inner  iteration  used  a  (preconditioned) 

m  — 

conjugate  gradient  method,  starting  with  the  zero  vector,  to  approxi- 
mately solve  a  system  of  the  form  Mz  =  r ,  in  the  sense  of  (1.2)  and 
(1.3).   The  symmetric  splitting  M  =  M,-M_  was  exploited  so  that  each 
step  of  the  inner  iteration  used  a  fast  direct  method  to  solve  a  system 
of  the  form  M^z"  =  r,  i.e.  the  Poisson  equation. 

In  addition  to  the  Richardson  results.  Table  1  also  gives  the 
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results  of  experiments  using  a  (preconditioned)  conjugate  gradient 
method  for  the  outer  iteration  to  solve  (3.2) .   See  Concus  and  Golub 
(1976)  and  Widlund  (1978)  for  a  discussion  of  this  method,  which  exploits 
the  symmetric/skew-symmetric  splitting.   The  inner  iterations  were 
carried  out  exactly  as  above,  being  terminated  by  (1.2)  and  (1.3) .   The 
conjugate  gradient  method  was  also  run  with  the  inner  iterations  carried 
out  to  machine  precision,  in  order  to  obtain  Oj  from  the  resulting  tri- 
diagonal  Lanczos  matrix.   These  runs,  made  solely  to  obtain  the  "optimal" 
CO  defined  by  (1.21),  are  not  included  in  the  table. 


15 


31 


Table  1. 

6_ 

METHOD 

1.   (RICHARDSON) 

m 

W 

001 

43 

332 

01 

44 

261 

1 

44 

169 

2 

47 

137 

5 

62 

117 

6 

71 

115 

7 

75 

109 

8 

147 

147 

001 

57 

446 

01 

46 

270 

1 

48 

180 

2 

52 

152 

5 

68 

132 

6 

78 

127 

7 

86 

132 

8 

162 

162 

CONJUGATE  GRADIENT 


m 

W 

36 

266 

38 

218 

52 

182 

— 

>500 

39 

42 
71 


279 
235 
237 
>500 


Note:  t  is  the  number  of  mesh  points  in  each  direction,  6  defines  (1.4)  , 
m  is  the  number  of  outer  iterates  and  W  is  the  total  number  of  inner 
iterates  (number  of  calls  to  the  fast  Poisson  solver) . 

The  results  of  Table  1  are  graphed  in  Figure  1.   Notice  that   the 
optimal  choice  of  6  ,  given  lo  and  a,  is  quite  well-determined  by  the 
curves,  and  has  a  value  which  can  be  close  to  one.   Using  a  near  optimal 
value  of  6  is  far  more  efficient  than  carrying  out  the  inner  iterations 
to  full  accuracy.   It  is  interesting  to  note  that  the  Richardson  method 
(with  near  optimal  oj  ,  in  practice  unknown)  is  less  efficient  than  the 
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Figure  1. 
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conjugate  gradient  method  when  the  inner  iterations  are  carried  out 
quite  accurately,  but  more  efficient  when  a  near  optimal  6  is  used. 

The  numerical  results  were  obtained  using  a  VAX-11/780  at  the 
Courant  Mathematics  and  Computing  Laboratory.   Double  precision  arith- 
metic, i.e.  approximately  16  decimal  digits  of  accuracy,  was  used. 
More  extensive  experiments  have  been  made,  but  they  will  be  reported 
in  a  subsequent  paper. 
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